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We report results of a numerical-relativity simulation for the merger of a black hole-neutron star 
binary with a variety of equations of state (EOSs) modeled by piecewise polytropes. We focus, in 
particular, on the dependence of the gravitational waveform at the merger stage on the EOSs. The 
initial conditions are computed in the moving-puncture framework, assuming that the black hole is 
nonspinning and the neutron star has an irrotational velocity field. For a small mass ratio of the 
binaries (e.g., Mbh/Mns = 2, where Mbh and Mns are the masses of the black hole and neutron 
star, respectively), the neutron star is tidally disrupted before it is swallowed by the black hole 
irrespective of the EOS. Especially for less-compact neutron stars, the tidal disruption occurs at a 
more distant orbit. The tidal disruption is reflected in a cutoff frequency of the gravitational- wave 
spectrum, above which the spectrum amplitude exponentially decreases. A clear relation is found 
between the cutoff frequency of the gravitational-wave spectrum and the compactness of the neutron 
star. This relation also depends weakly on the stiffness of the EOS in the core region of the neutron 
star, suggesting that not only the compactness but also the EOS at high density is reflected in 
gravitational waveforms. The mass of the disk formed after the merger shows a similar correlation 
with the EOS, whereas the spin of the remnant black hole depends primarily on the mass ratio of 
the binary, and only weakly on the EOS. Properties of the remnant disks are also analyzed. 

PACS numbers: 04.25.D-, 04.30.-w, 04.40.Dg 



I. INTRODUCTION 



Gravitational-wave observation is becoming one of 
the reliable tools for observing our Universe. Cur- 
rent ground-based laser-intcrfcromctric gravitational- 
wave detectors such as LIGO [l[ and VIRGO Q have 
already achieved nontrivial scientific results; e.g., upper 
limits on the amplitude of a stochastic gravitational- wave 
background have been improved and we now know that 
gravitational waves are not the main energy source of 
our Universe Q . Advanced gravitational- wave detectors 
such as advanced LIGO will be in operation within the 
next several years and detect gravitational waves, which 
can be used to explore the nature of strongly gravitat- 
ing phenomena. The most promising sources of gravi- 
tational waves are the coalescing compact binaries com- 
posed of compact objects such as black holes (BHs) and 
neutron stars (NSs). As illustrated in this paper, black 
hole-neutron star (BH-NS) binaries are potential sources 
for exploring the nature of the NSs and high-density nu- 
clear matter. 

According to a statistical study based on population 
synthesis calculations, the detection rate of gravitational 
waves from BH-NS binaries is estimated to be 0.5-50 
events per year for advanced gravitational-wave detec- 
tors 0, Q . This suggests that we will observe a variety of 
BH-NS binaries in the next decade. To extract physical 
information of BH-NS binaries as well as the informa- 
tion about the BH and NS themselves from gravitational 
waves, theoretical templates of gravitational waves are 



necessary. This fact motivates the numerical-relativity 
community to study in depth the coalescence of BH- 
NS binaries, because numerical relativity is the unique 
approach for accurately computing gravitational waves 
emitted from the late inspiral and merger phases of such 
compact binaries. 

Another astrophysical interest in BH-NS binaries is 
motivated by their potential to be a progenitor of short- 
hard gamma-ray bursts (GRBs); see and references 
therein for reviews. According to a merger scenario of 
GRBs, a NS is tidally disrupted by a low-mass BH be- 
fore the orbit reaches an innermost stable circular orbit 
(hereafter ISCO), resulting in a system consisting of a 
rotating BH and a hot, massive accretion disk of mass 
> O.O1M which could become the central engine of a 
GRB. This BH- massive disk system could subsequently 
radiate a large amount of energy > 10 48 erg by neutrino 
emission or by the so-called Blandford-Znajek process Q 
in a short time scale < 1 s. Then, neutrino-anti neutrino 
pair annihilation or electromagnetic Poynting flux could 
drive a GRB. One of the key questions for the merger 
scenario is whether the tidal disruption could lead to for- 
mation of the BH-disk system. Numerical relativity is 
again the unique approach for answering this question. 

The equation of state (EOS) of NSs, which is still un- 
known, is the key for determining gravitational wave- 
forms emitted in the tidal-disruption phase as well as for 
determining other properties of the BH-disk system such 
as the mass and typical density of the disk. The EOS 
(specifically its stiffness) determines the relation between 
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mass and radius of a NS, and hence, the relation between 
the tidal-disruption process and associated gravitational 
waveforms. The reason is that the sensitivity of a NS to 
the tidal force by the companion BH depends on its ra- 
dius; e.g., a NS of larger radius (with a stiffcr EOS) will 
be disrupted at a larger orbital separation (or a lower or- 
bital frequency). If the tidal disruption of a NS occurs at 
a larger distance, more material may be spread around 
the companion BH, and consequently, a high-mass rem- 
nant disk may be formed. Also, the gravitational-wave 
frequency at the tidal disruption, which will be one of 
the characteristic frequencies, is lower for NSs of larger 
radius. The EOS of nuclear matter beyond the normal 
nuclear density is highly uncertain due to the lack of con- 
straints obtained from experiments. Gravitational-wave 
astronomy will become a new and robust tool for deter- 
mining or at least constraining the EOS at such high 
densities through the observation of NSs [||-[l2j . For this 
purpose, we need theoretical templates of gravitational 
waves and it is necessary to perform many simulations 
employing a wide variety of possible EOSs for the NS 
matter. 

In recent years, fully general relativistic studies of BH- 
NS binaries have been perform ed both in calculations of 
quasiequilibrium states 



tions of the mergers [li 
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]] and in dynamical simula- 
However, we have not yet 



understood the effect of the EOS on the merger of BH-NS 
binaries in spite of its importance; in most of the previous 
studies, NSs are modeled by simple and unrealistic T = 2 
polytropic EOS (but see [25| )■ One of the next goals in 
numerical relativity is to clarify the effect of the EOS on 
the merger process of BH-NS binaries and on resulting 
gravitational waveforms. For this purpose, a systematic 
parametrization of possible EOSs by a small number of 
parameters is quite useful. 

In this paper, we report new results obtained by a sim- 
ulation using a wide variety of piecewise polytropic EOSs, 
which are shown to be useful for parametrizing nuclear- 
theory-based EOSs in the cold approximation [Ly, |26|, |27| 
[28}. We employ eight types of the piecewise polytropic 
EOSs, ranging from highly stiff to soft ones [29]. We 
systematically choose the BH and NS masses in a real- 
istic range of interest. As a first step in this series of 
work, the BH is assumed to be nonspinning. We track 
orbital evolutions of BH-NS binaries typically for ~ 5 
orbits so that the orbital eccentricity would not give a 
serious error in gravitational waveforms at the onset of 
the merger phase. We clarify the dependence of gravita- 
tional waveforms and merger remnants on the EOS. In 
particular, we show that a gravitational-wave spectrum 
contains valuable information on the EOS properties. 

This paper is organized as follows. In Sec. Ill Al we 
summarize initial conditions employed in this paper. Sec- 
tion flTB] describes the piecewise polytropic EOS and the 
models adopted in this paper. Section IIIII describes the 
formulation and methods of numerical simulations. Sec- 
tion IIVI presents the numerical results and clarifies the 
effect of the EOS on gravitational waveforms and merger 



remnants. Section IVl is devoted to a summary. Through- 
out this paper, we adopt the geometrical units in which 
G = c = 1, where G and c are the gravitational con- 
stant and the speed of light, respectively. The irreducible 
mass of the BH, gravitational mass of the NS in isolation, 
circumferential radius of the NS in isolation, Arnowitt- 
Deser-Misner (ADM) mass of the system, and sum of 
the BH and NS masses at infinite separation are denoted 
by A/ B h,Mns,-Rns,-M", m = A/ B h + M ns , respectively. 
The mass ratio Q is defined by Q = A/bh/A/ns an d the 
compactness of the NS (C) is defined by C = A/ns/^ns- 
Latin and Greek indices denote spatial and spacetime 
components, respectively. 



II. INITIAL CONDITION 

We employ BH-NS binaries in quasiequilibria for ini- 
tial conditions of numerical simulations as in [2(| l24j . 
The quasiequilibrium state is computed in the moving- 
puncture framework [l7l - [T9j with a piecewise polytropic 
EOS [ll|, [26|, [2?J ■ Here, we first summarize the formula- 
tion and numerical methods for computing the quasiequi- 
librium state and then describe EOSs employed in this 
paper. The details of the formulation and methods 
for computing initial conditions are described in [l7| . 
to which the reader may refer. Computation of the 
quasiequilibrium state is per formed using the spectral- 
method library LORENE [30j]. 



A. Formulation and methods 

We derive quasiequilibrium states of BH-NS binaries 
as solutions of the initial value problem of general rela- 
tivity [3lJ. When the orbital separation of the binary is 
large enough, the time scale for the gravitational-wave 
emission, £gw, is much longer than the orbital period 
P or b, so that we can safely neglect the radiation reaction 
of the gravitational-wave emission. In numerical simu- 
lations, the orbital evolution has to be followed for > 5 
orbits to derive a realistic waveform both for the late in- 
spiral and merger phases. For such a purpose, we have to 
choose the initial separation of the binary which satisfies 
*gw 3> -Porb, and have to provide BH-NS binaries in a 
quasicircular orbit as the initial condition, i.e., the binary 
is approximately in an equilibrium state if it is observed 
in the comoving frame. To satisfy these two conditions, 
we assume the presence of a helical Killing vector field 
with the orbital angular velocity fl, 



(i) 



and a hydrostatic equilibrium for the fluid configuration 
in the comoving frame. In addition, we assume that the 
BH is nonspinning and the NS has an irrotational ve- 
locity field. The irrotational velocity field is believed to 
be an astrophysically (approximately) realistic configu- 
ration [32, 1321 . 
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To compute the three- metric 7^, the extrinsic curva- 
ture Kij, the lapse function a, and the shift vector ft' 1 , 
we employ a mixture of the conformal thin-sandwich ap- 
proach and the conformal transverse-traceless decompo- 
sition of Einstein's equation [3l|. We assume the con- 
formal flatness of the three- metric 7y — ip 4 ^fij — ip 4 fij, 
the stationarity of the conformal three-metric dt'y-ij = 0, 
and the maximal slicing condition for the trace part of 
the extrinsic curvature K = jijK^ , i.e., K = d t K = 0. 
Here, /y denotes the flat spatial metric. Then, the basic 
equations for the conformal factor tp, the shift vector ft 1 , 
and a weighted lapse function $ = aip are derived from 
the Hamiltonian constraint, the momentum constraint, 
and the maximal slicing condition dtK = as 

Aib = -2n^ 5 p H - l^AijA^ , (2) 
8 

Aft 1 + \v l V 3 ft J = 167r$V 3 / + 2A tl V J (^~ 7 ), (3) 

A$ = 2tt$V z W + 25) + l^-^A^A^, (4) 

8 

where A ij = ib w K^, A = /^ViVj, and V* denotes the 
covariant derivative associated with /y. We assume an 
ideal fluid for the matter field 

= phu»u u + Pg^. (5) 

where p is the rest-mass density, P is the pressure, h = 
1+e+P/p is the specific enthalpy, e is the specific internal 
energy, and is the four- velocity of the fluid. Then, the 
fluid quantities seen by the Eulerian observer are denoted 
by 

p H = phiau 1 ) 2 - P, (6) 
f = phau'u^;, (7) 
S = ph^au 1 ) 2 - 1] + 3P. (8) 

The EOS fully determines relations among the thermo- 
dynamical quantities p, e, P, and h. We describe the 
EOS adopted in this work in Sec. Ill Bl 

In the moving-puncture framework, we set ip and <E> as 

^ = l + -^ + 0,$ = l--^ + ,, 0) 

2r B H ?"BH 

where Mp and M$ are positive constants of mass dimen- 
sion and tbh = \x* — x\,\ is a coordinate distance from 
the puncture located at x\>. We numerically solve the 
nonsingular parts <fi and rj using Eqs. @ and ^ and ad- 
justing the parameter Mp to achieve a desired BH mass. 
The other parameter, M$, is determined by the virial re- 
lation, i.e., the condition in which the ADM mass (Mo) 
and the Komar mass agree, which holds for the station- 
ary and asymptotically flat spacetime [U, [H| , 

<b di$dS l = - f d l ipdS i = 2ttM . (10) 

We note that the lapse function, a, obtained in this 
method is always negative near the puncture. In the 



numerical simulation, we modify the initial condition for 
a appropriately and ensure its positivity. 

For solving the momentum constraint, we decompose 

j j ctS 

kj = ViWj + VjWi - p tJ V k W k + Kfj, (11) 

where Wi is an auxiliary three-vector field, and W 1 = 
pi Wj . Kfj denotes a conformally weighted extrinsic cur- 
vature associated with the linear momentum of the BH, 
written by [36| 

K i - ^-[^ BH + hp? n (fij - hh)i k pri (12) 

z 'bh 

where P = £ bh /Vbh is a unit radial vector, 1^ = fijl 3 , 
and -Pj BH is the linear momentum of the BH, which is 
determined by the condition in which the total linear 
momentum of the system should vanish, 

ff H =- j 'j^ 6 d 3 x. (13) 

Wi obeys an elliptic equation 

AWi + ^ViVjVF = 87r^ 6 ii, (14) 

which is derived by taking a derivative of Eq. ([lip and 
using the momentum constraint. 

To summarize, we solve the elliptic equations for (j), ft 1 , 
77, and Wi imposing outer boundary conditions derived 
from the asymptotic flatness. In the present formalism, 
we do not have to impose inner boundary conditions at 
the BH horizon unlike in the excision method [HI, EE] • 

The basic equations for the hydrostatic equilibrium are 
derived from the condition of irrotation, i.e., the zero 
relativistic vorticity 

= 0, (15) 

and from the helical symmetric relation for the specific 
momentum of the fluid £^(hu^) = 0. One result is the 
first integral of the relativistic Euler equation, 

= -C(= const). (16) 

This equation determines h (and subsequently p, e, and 
P through an EOS) for an arbitrarily chosen constant C. 
The irrotational flow condition implies the presence of a 
velocity potential 'J, which determines the four- velocity 
of the fluid by hui — D^, where Di is the covariant 
derivative associated with 7^ . The continuity equation 
V Al (pu M ) = then leads to an elliptic equation for the 
velocity potential 

A quasiequilibrium state is computed using an itera- 
tion method described in detail in (l7j . During the it- 
eration, we fix the center of mass of the binary with a 
3PN- J condition described in [HI l2~i| ; we determine the 
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center of mass phenomenologically so that the total an- 
gular momentum of the binary for a given value of Ojtio 
agrees with that derived by the third post-Newtonian 
(3PN) approximation [37j ■ In this condition, the initial 
orbital eccentricity is by a factor of ~ 2 smaller than that 
in other conditions tried to this time [l8T[20l | , and the ec- 
centricity at the onset of the merger becomes < 1% for a 
longterm simulation which tracks ~ 5 inspiral orbits. 



B. Piecewise polytropic equation of state 

The temperature of NSs, except for newly born ones, 
are believed to be much lower than the Fermi energy of 
the constituent particles (38[- This implies that we can 
safely neglect the thermal effects and employ a cold EOS, 
for which the pressure, P, the specific internal energy, e, 
and other thermodynamical quantities are written as a 
function of the rest-mass density p. One of the simplest 
cold EOSs is a polytropic EOS, 



P = np 1 



-l/r. 



(17) 



where k is the polytropic constant and n p (> 0) the poly- 
tropic index: In the following, we often refer to the adia- 
batic index defined by T = 1 + l/n p . The first law of the 
thermodynamics, 



ds 



P<1 [ - 

,P 



(18) 



determines the specific internal energy as e — P/[(T— l)p] 
where we assume e = at p = 0. Then, the specific 
enthalpy h becomes 



h = 1 



r 



„r-i 



(19) 



A piecewise polytropic EOS is a phenomenologi- 
cally parametrized EOS, which reproduces cold nuclear- 
theory-based EOSs at high density only with a small 
number of polytropic constants and indices [HI HH, H3] , 
i.e., 

P(p) = K t p r - for < p < Pi (1 < i < n), (20) 

where n is the number of the pieces used to parametrize 
an EOS and pi denote boundary densities for which we 
provide an appropriate value (see the method below). 
Here, po = and p n — ¥ oo. Ki is the polytropic constant 
and Ti the adiabatic index for each piece. We note that 
we could in principle match to the known more realistic 
EOS at lower density. However, using a single polytrope 
for the low-density EOS is justified to the extent that the 
radius and deformability of the NS as well as resulting 
gravitational waveforms in the merger phase are insensi- 
tive to the low-density EOS. 

At each boundary density, p — pi (i = l,...,n— 1), 
the pressure is required to be continuous, i.e., Kip i i — 
Ki + ip\ %+1 . Thus, if we give K\, Ti, and pi (i = 1, n), 



the EOS is totally determined. For the zero-temperature 
EOS, the first law of the thermodynamics (fT5|) holds, 
and thus, e and h are also determined except for the 
choice of the integration constants, which are fixed by 
the continuity condition of e (hence equivalently h) at 
each pi. 

Recently, several authors have shown that the piece- 
wise polytropic EOS composed of one piece in the crust 
region and three pieces in the core region approximately 
reproduces most of nuclear-theory-based EOSs at high 
density [26| . Here, three pieces in the core region are re- 
quired to reproduce a high-mass NS for which inner and 
outer cores could have different stiffness due to the vari- 
ation of properties of high-density nuclear matter. In the 
present work, we pick up NSs of relatively low mass 1.2- 
1.35Af Q , taking into account that the masses of the NSs 
in the observed binary are fairly small (39|. The highest 
density of such NSs is not high enough in general that 
the EOS for the high-density part plays a critical role 
(note that if the EOS is very soft, the EOS for the high- 
density region is important, but we do not pursue this 
possibility in this paper) . An additional fact to be noted 
is that NSs in BH-NS binaries never achieve the state of 
density higher than the initial value; their density should 
decrease due to the tidal field of the companion BH dur- 
ing the evolution. For these reasons, we employ a sim- 
ple version of piecewise polytropic EOS in this paper, in 
which only one piece is assigned for the core region and 
one piece for the crust region as in [111 ], Following [Tl| . 
we employ the parameters of the crust EOS for all the 
models as follows: 



r a = 1.35692395, 

ki/c 2 = 3.99873692 x 1CT 8 £ 



i-r lcm3ri 



(21) 
(22) 



On the other hand, we vary the value of T 2 (the adiabatic 
index for the core EOS). Authors in [ll[ propose that in- 
stead of giving the density p\ , the pressure p at the fidu- 
cial density p^du = 10 14 7 g/cm 3 in the core region should 
be provided because this parameter p is closely correlated 
with the NS radius and deformability [40] . Thus, we have 
the following relations: 



P - 



(= P(Pl))- 



(23) 
(24) 



These determine the values of k 2 and p\. 

Table Q] lists the parameters of the EOSs employed 
in this paper, and several key quantities for each EOS. 
"2H," "H," "HB," and "B" denote very stiff, stiff, moder- 
ately stiff, and soft EOSs, respectively, for which T 2 = 3.0 



universally, but the values of p are varied [llj. For "HB," 



"HBs," and "HBss" or "B," "Bs," and "Bss", we assign 
the same value of p but different values of T 2 . The sub- 
script "s" denotes that the value of T 2 is smaller. For "s" 
and "ss," T 2 = 2.7 and 2.4, respectively. 

We calculate all the physical quantities for the spher- 
ical NS in equilibrium both by solving the Tolman- 
Oppenheimer-Volkoff equation directly and using the 
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code to calculate initial conditions in the isotropic gauge 
by LORENE, and check that numerical values agree with 
each other within 0.03%. Figure Q] plots the relation be- 
tween the mass Mns and circumferential radius -Rns for 
the spherical NSs with the adopted piecewise polytropic 
EOSs. For comparison, we also plot the relation for T = 2 
polytropic EOS with k/c 2 = 2 x 10 _16 g _1 cm 3 . Note 
that in the polytropic EOS with a fixed adiabatic index, 
only the shape of this relation has an invariant meaning 
and there is a freedom of the absolute scaling, since all 
the dimensional quantities can be rescaled through the 
polytropic length scale i? po i y = K 1 '( 2r ~ 2 '. 

Figure [T] shows that for a given mass ~ 1.35M0, the 
radius depends strongly on the EOSs, whereas the radius 
for a given piecewise polytropic EOS depends only weakly 
on the mass around the canonical mass ~ 1.35M©. This 
weak dependence of the radius on the mass is an often- 
seen feature for the nuclear-theory-based EOSs (38|. By 
contrast, the relation calculated with the T = 2 poly- 
tropic EOS does not show this feature. Figure [T] illus- 
trates that the dependence of the radius i?NS on the 
mass Mns becomes much stronger in this EOS than in 
the piecewise polytropic EOSs. This illustrates that the 
r = 2 polytropic EOS is not very realistic. 

Comparison of the quantities among HB, HBs, and 
HBss EOS models in Table Q] reveals a complicated mass- 
radius relation: HB is not always stiffer than HBss. In- 
deed, the radius with Mns = 1-2M© is largest for HBss 
and smallest for HB among three models, whereas the ra- 
dius with Mns = 1.35M is largest for HB and smallest 
for HBss. This complicated relation of the "stiffness" is 
due to the choice for the combination , p) (cf. Tablefl]) . 
For a density smaller than p^du, HBss EOS is stiffer than 
HB and HBs EOSs, whereas for a high density p > padu, 
HB EOS is stiffer than the others. For a given high-mass 
NS for which the central density is much larger than pfidu, 
the radius with HB EOS should be larger than that with 
other two EOSs. By contrast, for a given low-mass NS 
for which the central density is not very high, the radius 
with HB EOS should be smallest. 



C. Models 

The previous works by three groups [22l42^ | have found 
that the NSs in BH-NS binaries with high mass ratio 
Q > 4 are barely subject to tidal disruption if the com- 
panion BH is not spinning: At the merger, the BH swal- 
lows most of the NS matter at one moment and the rem- 
nant disk mass is quite small or nearly equal to zero. 
Namely, the NS behaves approximately as a point par- 
ticle even at the ISCO. Gravitational waves emitted in 
such a case have a similar waveform to that from a BH- 
BH binary. Because the behavior of NSs with high-mass 
BH companions does not show remarkable dependence 
on the EOS, they are unsuitable for the purpose of this 
paper, i.e., to investigate the effect of the EOS on grav- 
itational waves and final outcomes. Thus, we focus only 
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FIG. 1. The relation between the mass and circumferential 
radius of the spherical NSs for piecewise polytropic EOSs 
adopted in this paper. For comparison, we also plot the curve 
for r = 2 polytropic EOS with s/c 2 =2x lGr 16 g _1 cni 3 (dot- 
ted curve). 

on low mass-ratio binaries with Q = 2 and 3 in this pa- 
per. Also, we choose relatively low-mass NSs, because 
two-piece EOSs adopted in this paper may not be ap- 
propriate for modeling a high-mass NS with high central 
density, due to the lack of model parameters in the high- 
density region. 

Table ITT1 summarizes key quantities for the initial mod- 
els employed in the present numerical simulation. The 
labels for the models denote the name of the EOS, the 
mass ratio, and the NS mass; e.g., 2H-Q2MI35 is mod- 
eled by 2H EOS, and its mass ratio and the NS mass are 
2 and 1.35M©, respectively. The primary purpose of this 
paper is to study the dependence of gravitational wave- 
forms and the final outcome on (i) the EOS of NSs, (ii) 
the mass ratio, and (hi) the NS mass. These purposes 
are reflected in our choice of the initial models. 

We prepare quasiequilibrium states basically with the 
same value of f2mo for the same value of Q irrespective of 
the EOS. The value of fimo is chosen to be small enough 
that the binaries spend more than 5 inspiral orbits before 
the onset of the merger. For Q = 2 binaries, a smaller 
value of initial angular velocity is required only for 2H 
EOS, because the NS with this EOS has a much larger 
radius than with other EOSs and is sensitive to the BH 
tidal force even for a larger orbital separation; to track 
> 5 inspiral orbits before the tidal disruption, we have 
to choose the value of Qm by ~ 10% as small as that 
for other EOSs. For the case of Q = 3, we also choose 
smaller values of fimo for Mns = 1.2M© cases. 



III. METHODS OF SIMULATIONS 

Numerical simulation is performed u sing an adaptive- 
mesh refinement (AMR) code SACRA [4jJ. The formu- 
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TABLE I. Key ingredients of the adopted EOSs. T2 is the adiabatic index in the core region and p is the pressure at the fiducial 
density padu = 10 14 ' 7 g/cm 3 , which determines the polytropic constant K2 of the core region and p\\ the critical rest-mass 
density separating the crust and core regions. M max is the maximum mass of the NS for a given EOS. -R135 (-R12) and C135 (C12) 
are the circumferential radius and the compactness of the NS with Mns = 1.35M (1.2Mq). 



Model 


r 2 


logioP (g/cm 3 ) 


pi (10 14 g/cm 3 ) 


M max [M ] 


-R135 (km) 


C135 


R12 (km) 


C12 


2H 


3.0 


13.95 


0.7033 


2.835 


15.23 


0.1309 


15.12 


0.1172 


H 


3.0 


13.55 


1.232 


2.249 


12.27 


0.1624 


12.25 


0.1447 


HB 


3.0 


13.45 


1.417 


2.122 


11.61 


0.1718 


11.60 


0.1527 


HBs 


2.7 


13.45 


1.069 


1.926 


11.57 


0.1723 


11.67 


0.1519 


HBss 


2.4 


13.45 


0.6854 


1.701 


11.45 


0.1741 


11.74 


0.1509 


B 


3.0 


13.35 


1.630 


2.003 


10.96 


0.1819 


10.98 


0.1614 


Bs 


2.7 


13.35 


1.269 


1.799 


10.74 


0.1856 


10.88 


0.1629 


Bss 


2.4 


13.35 


0.8547 


1.566 


10.27 


0.1940 


10.66 


0.1663 



TABLE II. Key parameters and quantities for the initial conditions adopted in the numerical simulations. The adopted EOS, 
mass ratio (Q), NS mass in isolation (Mns), angular velocity (fi) in units of c 3 /Gmo, baryon rest mass (M 4 ), compactness of 
the NS in isolation (C), maximum rest-mass density (p max ), ADM mass of the system (Mo), and total angular momentum of 
the system (Jo), respectively. 



Model 


EOS 


Q M NS [M e ] 


GQm /c 3 


M*[M Q ] 


C 


pmax ( 


g/cm 3 ) 


M [Mq\ 


MGM 2 @ /c] 


2H-Q2M135 


2H 


2 


1.35 


0.0250 


1.455 


0.1309 


3.740 


xlO 14 


4.015 


14.39 


H-Q2M135 


H 


2 


1.35 


0.0280 


1.484 


0.1624 


7.018 


xlO 14 


4.013 


14.02 


HB-Q2M135 


HB 


2 


1.35 


0.0280 


1.493 


0.1718 


8.262 


xlO 14 


4.013 


14.02 


HBs-Q2M135 


HBs 


2 


1.35 


0.0280 


1.489 


0.1723 


9.154 


xlO 14 


4.013 


14.02 


HBss-Q2M135 


HBss 


2 


1.35 


0.0280 


1.485 


0.1741 


1.082 


xlO 15 


4.013 


14.02 


B-Q2M135 


B 


2 


1.35 


0.0280 


1.503 


0.1819 


9.761 


xlO 14 


4.013 


14.02 


Bs-Q2M135 


Bs 


2 


1.35 


0.0280 


1.501 


0.1856 


1.137 


xlO 15 


4.013 


14.02 


Bss-Q2M135 


Bss 


2 


1.35 


0.0280 


1.501 


0.1940 


1.490 


xlO 15 


4.013 


14.02 


2H-Q3M135 


2H 


3 


1.35 


0.0280 


1.455 


0.1309 


3.737 


xlO 14 


5.359 


21.05 


H-Q3M135 


H 


3 


1.35 


0.0300 


1.484 


0.1624 


7.011 


xlO 14 


5.358 


20.74 


HB-Q3M135 


HB 


3 


1.35 


0.0300 


1.493 


0.1718 


8.254 


xlO 14 


5.358 


20.74 


B-Q3M135 


B 


3 


1.35 


0.0300 


1.503 


0.1819 


9.751 


xlO 14 


5.357 


20.74 


2H-Q2M12 


2H 


2 


1.20 


0.0220 


1.282 


0.1172 


3.466 


xlO 14 


3.571 


11.71 


H-Q2M12 


H 


2 


1.20 


0.0280 


1.303 


0.1447 


6.421 


xlO 14 


3.567 


11.08 


HB-Q2M12 


HB 


2 


1.20 


0.0280 


1.310 


0.1527 


7.522 


xlO 14 


3.567 


11.08 


B-Q2M12 


B 


2 


1.20 


0.0280 


1.317 


0.1614 


8.832 


xlO 14 


3.567 


11.08 


HB-Q3M12 


HB 


3 


1.20 


0.0280 


1.310 


0.1527 


7.517 


xlO 14 


4.763 


1.663 


B-Q3M12 


B 


3 


1.20 


0.0280 


1.317 


0.1614 


8.826 


xlO 14 


4.763 


1.663 



lation, the gauge conditions, the numerical scheme, and 
the methods of diagnostics are essentially the same as 
those described in [24| El| except for the EOS. Thus, 
we here only briefly review them. We also describe the 
present setup of the computational domain for the AMR 
algorithm and grid resolution in Sec. MI CI 



A. Formulation and numerical methods 

In SACRA, we solve Einstein's evolution equation in 
the BSSN formalism [42|, |43[ with the moving-puncture 
method [3(1 0, [45| . We evolve a conformal factor W = 
the conformal three-metric 7^ = 7 _1 ^ 3 7y, the 



trace of the extrinsic curvature K, the conformal trace- 
free part of the extrinsic curvature Aij = 7~ 1 / 3 (i"Q. ( — 
K-fij/3), and an auxiliary variable T l = —0^7^ . The spa- 
tial derivatives in the evolution equations are evaluated 
by a fourth-order centered finite difference except for the 
advection terms which is evaluated by a fourth-order non- 
centered finite difference. A fourth-order Runge-Kutta 
method is employed for the time evolution. 

Following [46j , we employ a moving-puncture gauge in 
the form 

(d t - P J d 3 )a = -2aK, (25) 
(8 t - pd 3 )P = (3/4)5\ (26) 

(d t - vih ur = (d t - ^a,)f* - Vs b\ (27) 
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where B l is an auxiliary variable and rj s is an arbitrary 
constant. In this work, we typically set r) s ps Mq/Mbh- 

For the hydrodynamics, we evolve p„ = pai^W' 3 , 
&i = /lUj, and e* = /iau* — P/(pau t ). To handle the ad- 
vection terms, we adopt a high-resolution central scheme 
by Kurganov and Tadmor [471 ] with a third-order piece- 
wise parabolic interpolation for the cell reconstruction. 

With regards to the EOS, we decompose the pressure 
and the specific internal energy into cold and thermal 
parts as follows (e.g., (48[) 



P — -Pcold + Pth , £ — Scold + £th- 



(28) 



Here, the thermal part is nonzero only in the presence 
of shock heating, and thus, this part plays a role for the 
evolution only in the merger phase. Once the primitive 
variables p and e are recovered from the conserved vari- 
ables p*, Mi, and e*, we calculate zero-temperature parts 
P C oid and e co id from p using the piecewise polytropic EOS 
(I2U1) . Then, the thermal part of the specific internal en- 
ergy is calculated by e t h = £ — e C oid: and finally the ther- 
mal part of the pressure Pth is determined. In this paper, 
we adopt a simple T-law, ideal-gas EOS for the thermal 
part as (e.g., |48j |) 



Pth = (I\h - l)p£th, 



(29) 



where I\h is an adiabatic index for the thermal part. 
We choose I\h equal to the adiabatic index in the crust 
region, Tx, for simplicity. 

Because the vacuum is not allowed in any conservative 
hydrodynamic scheme, an artificial atmosphere of small 
density is distributed outside the NS in the same manner 
as done in our previous work [24j]. The rest-mass den- 
sity of the atmosphere is set to be p atm — lO^Vmax ~ 
10 6 g/cm 3 for the inner computational domain. For the 
outer domain with r > r c w 20Pns, a smaller density is 
assigned according to the rule p — p a tme- 1 ~ r ^ rc ■ The total 
rest mass of the atmosphere is always less than 1O _5 M0, 
and hence, we can safely neglect spurious effects by accre- 
tion of the atmosphere onto the remnant accretion disk 
as long as the disk mass is much larger than 10 _5 M©. 



B. Diagnostics 

Gravitational waves are extracted calculating the out- 
going part of the complex Weyl scalar ^4, which we eval- 



uate at a finite coordinate radii r = 300-400M Q . Gravi- 
tational waveforms are obtained by integrating ^4 twice 
in time as 



h + (t)-ih x (t) 



dt' / dt"y A (t"), 



(30) 



and then by subtracting the quadratic function c^i 2 + 
a\t + ao from the obtained waveform using the least- 
square fitting for determining the constants ao, ax, and 
a-i- The purpose of this subtraction is to eliminate 
unphysical components in numerically calculated Weyl 
scalar, v& 4 [Hj], as described in [24| (see also [I(|). We 
also calculate the amount of radiated energy AE and 
angular momentum A J by integrating the emission rate 
calculated from the Weyl scalar ^4 as 



dE 
~dt 

dJ z 
dt 



r 

I671- j s 

—Re 
16tt 



^ 4 dt 



dA, 
^ 4 dt 

J J d^idtdt'j dA 



(31) 



, (32) 



where S denotes a coordinate sphere of r — const, 
dA = r 2 d(cos6)dip is the surface element of S, and ^4 
is the complex conjugate of ^4. We decompose ^4 into 
,s = —2 spin-weighted spherical harmonics of 2 < I < 4. 
Among them, (Z, |m|) = (2,2) modes are always domi- 
nant but higher I modes such as (I, |m|) = (3,3), (4,4) 
and (2, 1) modes contribute to the totally radiated energy 
and angular momentum by larger than 1%. 

We compare numerical waveforms with those derived 
by the Taylor-T4 formula in the post-Newtonian approx- 
imation [51j for two point masses in quasicircular orbits. 
Assuming that both the BH and NS have no spin angu- 
lar momentum, we calculate the evolution of the orbital 
angular velocity f2(i) through X(t) — [mo^(t)] 2 ^ 3 and 
the orbital phase 8(i) up to 3.5PN order by solving the 
ordinary differential equations 52] 



dX 
~~dt 



5mo 

/4159 

~ ^"672" 

541 2 

H v ■ 

896 



1 



743 + 924t/ 

336 
15876 



X + 4ttX 3/2 



/ 34103 13661 



672 
5605 



V 18144 2016 
16447322263 1712 



139708800 

/4415 



59 

I 1 

18 

16 2 

Tor* + 



56198689 451 



48 



217728 



2592 



105 v ') V 4032 



358675 
6048 



91495 
1512 1 



(33) 
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where v = Q/(l + Q) 2 and je is the Euler constant. 



After X and are obtained, we calculate complex 
gravitational- wave amplitude h 22 of (I, m) — (2, 2) mode, 
assuming that the binary is orbiting on the equatorial 
(9 = 7r/2) plane, up to 3PN order using the formula (53j 



wing - 2 »e Y 
5 L> 



34 



107 
~42~ 



55 i v 
42^ 1X 



^ 7T 



107 

~2T 2! 

278185 41 



33264 96 



24ii/| X 5 / 2 
20261 9 



2ttX 

J 27027409 
\ 646800 
114635 , 



3/2 



2173 1069 



\1512 

856 2 2 
105 7 * + 3* 



216 



2047 
1512 1 



1712, „ 428, v 

In 2 lnX 

105 105 



2772 



99792 



428 
105' 



in \X 



(35) 



where D is a distance between the binary and an ob- 
server. We also compute a gravitational-wave spectrum 
from the waveform obtained in this way. 

We determine the properties of the BHs formed after 
the merger such as masses and spins using the quanti- 
ties associated with the apparent horizon. The apparent 
horizon is determined in the same manner as described 
in [U. 

The BH mass may be estimated by two methods. In 
the first method, we measure the circumferential radius 
C e of the apparent horizon along the equatorial plane and 
calculate C e /4w, which gives the BH mass in the station- 
ary vacuum BH spacetime. By this method, we estimate 
the mass of the remnant BH after the spacetime settles to 
an approximately steady state (assuming that deviation 
from Kerr spacetime due to the presence of surround- 
ing materials is negligible). In the second method, we 
measure the irreducible mass of the BH, Mj rr , which is 
determined from the area of the apparent horizon Aah 
as 



M, 



A 



AH 



16tt 



(36) 



For the Kerr spacetime, M irr is written by the mass and 
dimensionless spin parameter a = Jbh/^bh °f a BH 
(where Jbh is the spin angular momentum of the BH) as 



Mw r = Ml 



BH* 



(37) 



Thus, if either the BH spin or mass is known, the BH 
mass or spin is determined (again assuming that devia- 
tion from Kerr spacetime due to the presence of surround- 
ing materials is negligible). For the Kerr spacetime, this 
relation may be written as 



M irr = 



C' 



(38) 



Thus, we may say that the spin is estimated by calculat- 
ing Mjrr and C e . 



The dimensionless spin parameter of the BH is esti- 
mated also using the quantities defined on the apparent 
horizon. For a Kerr BH with spin parameter a, the ratio 
of the circumferential radius along the meridional plane 
C p to the one along the equatorial plane C e is written as 



-E 



2f, 



(39) 



where f + = 1 + yl — a 2 is the normalized radius of the 
horizon and E{z) is an elliptic integral 



r /2 i 

E{z) = I VI -zsin 2 6d6. 
Jo 



(40) 



Assuming that this relation holds for a BH surrounded 
by materials again, we estimate the spin parameter of the 
remnant BH. 

Comparison of the spin obtained from C p /C e with that 
derived from Eq. (|38t provides a consistency check. It is 
found that these two values agree with each other within 
the error Aa = 0.003 irrespective of the model of BH- 
NS binaries. For this reason, in the following, we only 
present the spin determined from C p /C e . 

In addition to the quantities for the remnant BHs, we 
calculate the total rest mass of materials located outside 
the apparent horizon by integrating the rest-mass density 
with respect to the proper volume element, 



M 



r>r A H 



p*d x, 



(41) 



r>rAH 



where tah = ''ah {6, cp) denotes the radius of the apparent 
horizon as a function of the angular coordinates (9,ip). 
M r > rAH is regarded as the mass of the remnant disk when 
the system settles to a quasistationary state after the 
merger. 



C. Setup of AMR grids 

Numerical simulation is performed using an AMR al- 
gorithm described in |4lj |. to which the reader may refer 
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TABLE III. Setup of the grid structure for the computation 
with our AMR algorithm. Ax = h 6 = L/(2 e N) is the grid 
spacing at the finest-resolution domain with L being the loca- 
tion of the outer boundaries for each axis. i?diam/Ax denotes 
the grid number assigned inside the semimajor diameter of 
the NS. Ao is the gravitational wavelength of the initial con- 
figuration. 



ivioaei 


IAX j IVlQ 


7? / A t> 
-Kdiam/iiiE 


J-i/ Ao 


2H-Q2M135 


0.0471 


90.8 


1.189 


H-Q2M135 


0.0377 


86.2 


1.065 


HB-Q2M135 


0.0347 


87.0 


0.982 


HBs-Q2M135 


0.0353 


85.2 


0.998 


HBss-Q2M135 


0.0353 


84.0 


0.998 


B-Q2M135 


0.0330 


85.1 


0.932 


Bs-Q2M135 


0.0324 


84.4 


0.915 


Bss-Q2M135 


0.0270 


95.4 


0.825 


2H-Q3M135 


0.0353 


89.0 


0.998 


H-Q3M135 


0.0282 


84.7 


0.856 


HB-Q3M135 


0.0269 


82.7 


0.816 


B-Q3M135 


0.0247 


83.8 


0.749 


2H-Q2M12 


0.0565 


86.9 


1.255 


H-Q2M12 


0.0453 


83.1 


1.281 


HB-Q2M12 


0.0420 


83.6 


1.188 


B-Q2M12 


0.0392 


83.4 


1.109 


HB-Q3M12 


0.0306 


84.6 


0.866 


B-Q3M12 


0.0278 


86.9 


0.786 



for details. In the present work, we prepare seven re- 
finement levels to ensure that the computational domain 
extends to the local wave zone for initial quasiequilib- 
rium states and that both compact objects are resolved 
with a sufficient grid resolution (e.g., Table ITTT)) . Each re- 
finement domain consists of the uniform, vertex-centered 
grids with (2N + 1, 2N + l,N+l) grid points for (x, y, z) 
with the equatorial plane symmetry at z — imposed. In 
the present work, we typically choose N = 50, with the 
exception that N = 54 for model Bss-Q2M135, in which 
the NS is quite compact and needs to be resolved with 
a better grid resolution. For several models arbitrarily 
chosen, we performed numerical simulations with lower 
grid resolutions, N = 36 and 42, to check the conver- 
gence of the numerical results (see the Appendix). The 
edge length of the largest domain is denoted by 2L and 
the grid spacing for each domain is then hi — L/(2 l N), 
where I = 0-6. In all the simulations, two sets of four 
finer domains comoving with compact objects cover the 
region in the vicinity of two objects, and the other three 
coarser domains cover both objects by a wider domain 
with their origins being fixed at the approximate center 
of mass of the binary. Namely, we prepare 11 refinement 
domains in total for all the simulations. 

Table IIIII summarizes the parameters of the grid struc- 
ture for the simulations in this paper. As mentioned 
above, the value of L is chosen to be « Ao, where 




-60 -40 -20 20 40 60 
x [km] 

FIG. 2. Evolution of the coordinate separation of the binary 
a4p for model HB-Q2M135. 

Ao = tt/SIo is the gravitational wavelength at t = and 
f^o is the orbital angular velocity of the initial configu- 
ration. Because the gravitational wavelength decreases 
during the evolution of the binaries, the outer boundary 
of the computational domains is guaranteed to be located 
in the wave zone throughout the simulation. Each of the 
two finest domains covers the semimajor axis of the NS 
with 42-48 grid points and the BH radius (the coordinate 
radius of the apparent horizon) with typically ~ 20 grid 
points, respectively. For N = 54 run, the total memory 
required for the simulations is about 11.6 G bytes. We 
perform numerical simulations with personal computers 
of 12 G bytes memory and of core-i7X processors with 
clock speed 3.2 or 3.33 GHz. We only use two processors 
to perform one job with an OPEN-MP library. Typical 
computational time required to perform one simulation 
(for ~ 40 ms in physical time of coalescence) is 7-10 
weeks. 



IV. NUMERICAL RESULTS 

A. Orbital evolution and general merger process 

To obtain a realistic numerical result for gravita- 
tional waveforms and the final outcome formed after the 
merger, it is necessary to exclude spurious effects associ- 
ated with a noncircularity in the orbital motion as much 
as possible. To assess the circularity of the orbital mo- 
tion, we plot the evolution of the coordinate separation 
a4 p = x l NS ~x l BH for model HB-Q2M135 in Fig.[2j Here, 
the position of the maximum rest-mass density is identi- 
fied as the coordinate of the NS, x^ s , and the location 
of the puncture, x P , is the coordinate of the BH, x B}i . 
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FIG. 3. Time evolution of the orbital angular velocity O(t)mo 
for models 2H-Q2M135, H-Q2M135, HB-Q2M135, and B- 
Q2M135 as a function of a retarded time defined by Eq. (|43[) 
with an appropriate time shift. The dotted curve denotes 
the evolution of the orbital angular velocity calculated by the 
Taylor- T4 formula. 



This figure suggests that the orbital eccentricity appears 
to be low throughout the whole evolution. Because > 5 
orbits are tracked, the eccentricity, which is likely to be 
nonzero initially, should be suppressed by gravitational 
radiation reaction. We note that for all the models, sim- 
ilar trajectories are found. 

The coordinate separation shown above is a gauge- 
dependent quantity. To show a stronger evidence that 
the eccentricity is suppressed to a small level, it is better 
to plot a gauge- independent quantity. Figure |3] plots the 
evolution of the orbital angular velocity defined from the 
(I, m) = (2, 2) mode of * 4 by 



n(t) 1 M = ™ = 2)| 



2 | J # 4 (Z = m = 2)dt\ ' 



(42) 



for models 2H-Q2M135, H-Q2M135, HB-Q2M135, and 
B-Q2M135. Here, the horizontal axis is chosen to be an 
approximate retarded time defined by 



< rct =t-D- 2M Q \n(D/M ). 



(43) 



We here do not plot the curve after the onset of tidal 
disruption. For comparison, the angular velocity derived 
from the Taylor-T4 formula is also plotted. To align the 
curve in the inspiral phase for Q(i)mo < 0.05, we ap- 
propriately shift the time for each model. For t rot < 
ms, an unphysical (a junk wave) component contained in 
the initial data dominates the waveform, and hence, fi(i) 
derived from Eq. (|42[) does not give the angular velocity. 

Figure [3] shows that the angular velocity obtained in 
numerical simulations agrees with that by the Taylor- T4 
formula within a small modulation of Afi/O < 5% irre- 
spective of the models. With the fact that the orbital 
eccentricity is approximately estimated as e ~ 2A£!/3£! 



for e <C 1 , we conclude that the orbital eccentricity is sup- 
pressed within ~ 3%. Figure[3]also shows that the devia- 
tion from the Taylor-T4 result becomes remarkable in an 
earlier time for models with stiffer EOSs such as 2H and 
H EOSs. This is due to the fact that the tidal elongation 
and disruption of the NS occur at slightly earlier stages 
of the inspiral orbits for models with the stiffer EOSs. 
This illustrates the fact that the stiffness of the EOS is 
reflected clearly in the gravitational- wave frequency (and 
gravitational-wave phase) as a function of time. 

Figures 2] and [5] plot the snapshots of the rest-mass 
density profiles and the location of the apparent horizon 
on the equatorial plane at selected time slices for models 
2H-Q2M12 and B-Q3M135. Figure H illustrates the pro- 
cess in which the NS is tidally disrupted to form a disk 
surrounding the companion BH. In this case, the NS is 
disrupted far outside the ISCO and then forms a one- 
armed spiral arm with large angular momentum. As a 
consequence of the angular momentum transport in the 
arm, a large amount of materials spread outward and 
then form a disk around the BH. We will report more 
details about the remnant disk in Sec. lIVDl Figure [S] il- 
lustrates the case in which the NS is not tidally disrupted 
before it is swallowed by the BH. In this case, mass of 
the disk formed after the onset of the merger is negligibly 
small. 



B. Gravitational waveforms 

Figures [5] and [7] plot the (l,m) — (2,2), plus-mode 
gravitational waveforms obtained numerically (hereafter 
referred to as h + ). All the waveforms are shown for an 
observer located along the z axis (axis perpendicular to 
the orbital plane) and plotted as a function of a retarded 
time t re t. We plot the amplitude in a normalized form, 
Dh+frnQ, and the physical amplitude observed by an ob- 
server located at a hypothetical distance D = 100 Mpc. 

To validate the numerical waveforms, we compare 
them with the Taylor-T4 waveform, which is accurate 
up to 3.5PN order in phase and 3PN order in amplitude, 
with an appropriate time shift; the time shift is carried 
out to align the curve of Q(t) as performed in Sec. IIV Al 
Figures [5] and [7J show that these two waveforms agree 
with each other irrespective of models during the inspi- 
ral phase, except for 2-3 initial cycles. The reasons for 
this initial disagreement are that an approaching veloc- 
ity associated with gravitational radiation reaction is not 
taken into account in the initial data and also the initial 
condition does not exactly model a quasicircular state, 
because we do not fully solve Einstein's equation for de- 
riving it. 

The numerical waveforms in the merger phase also (but 
due to a physical reason) deviate from the Taylor-T4 
ones both in phase and amplitude, in particular for mod- 
els with stiff EOSs, e.g., 2H-Q2M135 and 2H-Q2M12. 
For such models, ringdown waveforms associated with 
the BH quasinormal mode are not seen in the merger 
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FIG. 4. Evolution of the rest-mass density profile in units of g/cm and the location of the apparent horizon on the equatorial 
plane for model 2H-Q2M135. The filled circles denote the regions inside the apparent horizons. The color panels on the 
right-hand side of each figure show log 10 (p). 



and ringdown phases, and instead, the gravitational- wave 
amplitude damps suddenly in the middle of the inspiral 
phase. The reason for this quick damping is that the NS 
is tidally disrupted by the companion BH at an orbit in 
the inspiral phase within one orbital period, and then, 
the disrupted material forms a relatively low-density and 
nearly axisymmetric matter distribution around the BH, 
suppressing time variation of a mass quadrupole moment. 
Because the gravitational- wave emission stops in the mid- 
dle of the inspiral motion, the maximum amplitude of 
gravitational waves is smaller for such a binary than for 
a binary with no tidal disruption, as shown in Fig. [5] 
All these facts illustrate that the finite size effect of the 
NS significantly modifies gravitational waves derived in 
the point-particle approximation (in the Taylor-T4 for- 
mula). On the other hand, ringdown gravitational waves 
are clearly seen for models with soft EOSs (for which tidal 
disruption does not occur) such as model B-Q3M135, in 
which the numerical and the Taylor- T4 waveforms are in 
more excellent agreement even in the late inspiral phase. 

Table HVl presents total radiated energy AE and angu- 
lar momentum A J carried away by gravitational waves. 
The contribution from all the I — 2-4 modes is taken into 
account for AE and A J. We estimate systematic errors 
in the presented values to be less than 10%, which are as- 
sociated mainly with the finite grid resolution and partly 



with the finite extraction radii (cf. the Appendix). We 
note that the (I, \m\) — (2,2) modes always contribute 
by > 90% to both for AE and A J. The fraction of these 
modes is larger for binaries composed of less-compact 
NSs, because only binaries which escape the tidal disrup- 
tion in the late inspiral phase can efficiently emit higher 
Z-mode gravitational waves. Among other modes, (3, 3) 
and (4, 4) modes constitute most of the remaining part 
of A J, whereas the order of magnitude of the (2, 1) mode 
is as large as that of the (4, 4) mode for AE. 

The numerical results shown in Table IIVI illustrate a 
quantitative dependence of gravitational-wave emission 
on the compactness of the NS: For a given mass ratio, 
gravitational-wave emission continues for a longer dura- 
tion and consequently total radiated energy and angular 
momentum are larger for binaries composed of more com- 
pact NSs. Comparison among the models with Q = 2 
and Mns = 1.35Af© and with the same initial value of 
Qmo shows that both AE/Mq and AJ / Jq are monotoni- 
cally increasing functions of the NS compactness C. This 
point is also recognized from Figs. [5] and [71 e.g., from 
the comparison among gravitational waves for models 
H-Q2M135, HB-Q2M135, and B-Q2M135 (note that for 
model 2H-Q2M135 the simulation is started from a lower 
value of Qmo and it is not suitable for this comparison). 
Table [IV] also shows that AJ/AE decreases as the EOS 
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FIG. 5. The same as Fig. g] but for model B-Q3M135. 



softens. This is due to the fact that AJ/AE sa m/fl for 
a given angular harmonic of m, and for a soft EOS, more 
radiation is emitted at large angular velocity, f2. 



C. Gravitational- wave spectrum 

Characteristic features of a gravitational waveform, 
such as characteristic frequencies and their dependence 
on the EOS, are well reflected in the Fourier spectrum. 
Figures [8l[T0l display gravitational-wave spectra for all 
the models with the mass ratio Q = 2 and the models 
with the mass ratio Q = 3 and the NS mass Mns = 
1.35M . As before [24j, we define the Fourier spectrum 
as a sum of each Fourier component of two independent 
polarizations of the (7, m) = (2, 2) mode as 




M/)| 2 + |M/)| 2 



27Ilft h A (t)dt, 



(44) 
(45) 



where A denotes two polarization modes, + or x . In cal- 
culating h(f) from a numerically obtained Weyl scalar, 
^4, we always omit the unphysical radiation compo- 
nent extracted at t lct < ms using a step function 
of retarded time as the window function so that the 



TABLE IV. Total radiated energy AE and angular momen- 
tum AJ carried away by gravitational waves. AE and AJ 
are normalized with respect to the initial ADM mass Mq and 
angular momentum Jo, respectively. We also show the ratio 
between A J and AE. 



Model 


AE/M {%) 


AJ/J (%) (AJ/J )/(AE/M ) 


2H-Q2M135 


0.55 


14 


26 


H-Q2M135 


1.1 


20 


18 


HB-Q2M135 


1.4 


22 


16 


HBs-Q2M135 


1.4 


22 


16 


HBss-Q2M135 


1.5 


23 


15 


B-Q2M135 


1.7 


24 


14 


Bs-Q2M135 


1.9 


25 


13 


Bss-Q2M135 


2.2 


27 


12 


2H-Q3M135 


0.64 


15 


23 


H-Q3M135 


1.4 


22 


16 


HB-Q3M135 


1.6 


23 


14 


B-Q3M135 


1.8 


24 


13 


2H-Q2M12 


0.41 


12 


30 


H-Q2M12 


0.74 


16 


21 


HB-Q2M12 


0.89 


18 


20 


B-Q2M12 


1.1 


20 


18 


HB-Q3M12 


1.2 
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B-Q3M12 
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FIG. 6. (l,m) = (2,2), plus-mode gravitational waveforms for models 2H-Q2M135, H-Q2M135, HB-Q2M135, HBs-Q2M135, 
HBss-Q2M135, B-Q2M135, Bs-Q2M135, and Bss-Q2M135. All the waveforms are shown for an observer located along the z 
axis (axis perpendicular to the orbital plane) and plotted as a function of a retarded time. For model 2H-Q2M135, the waveform 
is plotted as a function of t re t — 5 ms to align it with other waveforms (note that the initial value of Q, only for this model is 
smaller than those for other models). The left axis denotes the amplitude normalized by the distance from the binary D and 
the total mass mo. The right axis denotes the physical amplitude of gravitational waves observed at a hypothetical distance 
100 Mpc. The dotted curves denote the waveform calculated by the Taylor-T4 formula. 
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FIG. 7. The same as Fig. H but for models 2H-Q3M135, H-Q3M135, HB-Q3M135, B-Q3M135, 2H-Q2M12, H-Q2M12, HB- 
Q2M12, and B-Q2M12. Again, the waveform for model 2H-Q2M12 is plotted as a function of t re t — 9 ms. 
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FIG. 8. Spectra of gravitational waves from BH-NS binaries for Q — 2 and Mns = 1.35M0 with all the EOSs chosen in this 
paper. The bottom axis denotes the normalized dimensionless frequency fmo(= Gfmo/c 3 ) and the left axis the normalized 
amplitude fh(f)D/mo. The top axis denotes the physical frequency / in Hz and the right axis the effective amplitude fh(f) 
observed at a distance of 100 Mpc from the binaries. The short-dashed slope line plotted in the upper left region denotes a 
planned noise curve of the Advanced-LIGO [lj optimized for 1.4M0 NS-NS inspiral detection ("Standard"), the long-dashed 
slope line denotes a noise curve optimized for the burst detection ("Broadband"), and the dot-dashed slope line plotted in the 
lower right region denotes a planned noise curve of the Einstein Telescope ("ET") The upper transverse dashed line is the 
spectrum derived by the quadrupole formula and the lower one is the spectrum derived by the Taylor- T4 formula, respectively. 
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FIG. 9. The same as Fig. |8] but for Q = 2 and for 
Mns = 1.35M© and 1.2Mq. Only the normalized amplitude 
fh(f)D/mo as a function of the dimensionless frequency /mo 
is shown. 



spurious radiation component does not introduce un- 
physical oscillations in the gravitational-wave spectrum. 
The spectrum amplitude for a low-frequency region of 
/ w S7(£ ret = 0)/-7r changes slightly if we include the spu- 
rious radiation component. However, we believe that our 



use of the window function is physically reasonable [55[ . 
We always show the spectrum based on gravitational 
waves observed along the z axis (axis perpendicular to 
the orbital plane), which is the most optimistic direction 
for the gravitational- wave detection. (To obtain an aver- 
aged amplitude, we only need to multiply a factor of 0.4; 



e.g., see [52J.) Because the Fourier components of any 
dimensionless quantity have the dimension of time, we 
define a dimensionless effective amplitude fh(f). In the 
figure, we plot this quantity observed at a hypothetical 
distance 100 Mpc as a function of / (Hz) or a normalized 
amplitude fh(f)D/mo as a function of dimensionless fre- 
quency fm . 

Figure [5] plots gravitational-wave spectra for Q = 2 
and Mns = 1.35M with all the EOSs employed in this 
paper. For all these models, the total mass is univer- 
sally too = 4.05M Q , and thus, a nondimensional quan- 
tity, /too(= Gfmo/c 3 ), is plotted at the bottom and / in 
units of Hz is plotted at the top. Also, a normalized am- 
plitude, fh(f)D/mo, is plotted at the left side and fh(f) 
observed at a distance of 100 Mpc from the binary is at 
the right side. For comparison, we also plot the spectra 
derived from the quadrupole formula (e.g., (56|) and the 
Taylor-T4 formula (dashed curves). 

General qualitative features of the gravitational-wave 
spectrum by BH-NS binaries are summarized as follows. 
In the early stage of the inspiral phase, during which the 
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FIG. 10. The same as Fig. [8] but for Mns = 1.35Af@ and for Q = 2 and 3. The left panel shows the normalized amplitude 
fh(f)D/mo as a function of the dimensionless frequency /mo. The right panel shows the spectra observed at a distance of 100 
Mpc. The spectra derived from the quadrupole formula and the Taylor- T4 formula are plotted by the short-dashed (Q = 2) 
and long-dashed lines (Q = 3). 



orbital frequency is < 1 kHz and the PN point-particle 
approximation works well, the gravitational-wave spec- 
trum is approximately reproduced by the Taylor- T4 for- 
mula. For this phase, the spectrum amplitude of fh(f) 
decreases as f~ ni where n.; = 1/6 for / -C 1 kHz and 
the value of increases with / for / < 1 kHz. As 
the orbital separation decreases, both the nonlinear ef- 
fect of general relativity and the finite size effect of the 
NS come into play, and as a result, the PN point-particle 
approximation breaks down. If the tidal disruption sets 
in for a relatively large separation (e.g. for 2H EOS), 
the amplitude of the gravitational-wave spectra damps 
for a low frequency in the middle of the inspiral phase 
(before the ISCO is reached). By contrast, if the tidal 
disruption does not occur or occurs at a close orbit near 
the ISCO, the spectrum amplitude for a high frequency 
region (/ > 1 kHz) is larger than that predicted by the 
Taylor-T4 formula (i.e., the value of n; decreases). In 
this case, an inspiral-like motion continues even inside 
the ISCO for a dynamical time scale and gravitational 
waves with a high amplitude are emitted. As a result, 
fh(f) becomes a slowly varying function of / for 1 kHz 
S f i5 /cut, where / cut ~ 2-3 kHz is the so-called cut- 
off frequency which depends on the binary parameters as 
well as the EOS of the NSs. (A more strict definition of 
/ cu t will be given below.) A steep damping of the spec- 
tra for / > / cut is universally observed, and for softer 
EOSs with a smaller radius of NSs, the frequency of / cu t 
is higher. This cutoff frequency is determined by the 
frequency of gravitational waves emitted when the NS is 
tidally disrupted for the stiff EOSs or by the frequency of 
a quasinormal mode of the formed BH for the soft EOSs. 
Therefore, the cutoff frequency provides potential infor- 
mation for the EOS through the tidal-disruption event of 
the NSs, in particular for the stiff EOSs. 

Hereafter, we pay special attention to the cutoff fre- 



quency determined by the tidal disruption. It is natural 
to expect that the NS compactness C primarily deter- 
mines the cutoff frequency in the combination, / cu tmo, 
because the orbital angular velocity at the onset of mass 
shedding, i? s hod, is written as a function of Q and C as 

MM 

Qm oc C . (46) 

In fact, we found a qualitative correlation between C and 
/ cut mo in the previous work [24]. To reconfirm this, we 
first plot gravitational- wave spectra [fh(f)D/mo as a 
function of /mo] for Q = 2 with the different NS mass 
A/ns = 1.35M and 1.2M© in Fig.H This indeed shows 
/cut mo increases monotonically with C irrespective of the 
NS mass for the given mass ratio. 

Figure [TU] shows the gravitational- wave spectrum for 
M NS = 1.35M and for Q = 2 and 3. The left panel 
plots fh(f)D/mo as a function of /too and the right 
panel fh(f) as a function of / for D = 100 Mpc. This 
shows that dependence of / cu tmo on C for Q = 3 is 
weaker than for Q = 2. The reason for this is that 
the tidal effect is weaker for Q = 3, as discussed in 
Sec. IIV Dl (As later shown in Fig. [Til /cut for models 
H-Q3M135, HB-Q3M135, and B-Q3M135 are not deter- 
mined by the orbital frequency at tidal disruption but 
by the quasinormal-mode frequency of the remnant BH, 
which sets an approximate upper limit on the frequency 
of gravitational waves emitted in the merger.) Hence, 
the information of the EOS is not encoded in gravita- 
tional waves for Q = 3 as strongly as for Q = 2. The 
right panel shows that / cu t is between ~ 1 and 3 kHz 
depending weakly on the value of Q. 

To analyze the cutoff frequency quantitatively and to 
strictly study its dependence on EOSs, we perform a sys- 
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tematic fitting procedure. As in [24|, we fit all the spectra 
by a function with seven free parameters 

M/) = ^3PN(/)e- (/// '- rins 

+ ^ e -(///-)°- [i _ e-t///-)'-], (47) 

where /z.3pn(/) is the Fourier spectrum calculated by the 
Taylor- T4 formula and / ins , / ins2 , / cu t, 01ns, er ins2 , cr cut , 
and A are free parameters. The first and second terms of 
Eq. (|47|) denote the spectrum models for the inspiral and 
merger phases, respectively. We determine these free pa- 
rameters by searching the minimum for a weighted norm 
defined by 

£{[M(/0 /a. (/.)./•; 3 } 2 , (48) 

i 

where i denotes the data point for the spectrum. In the 
previous work (24j . we fix <7i ns = 3.5 and <7i nS 2 = 5 to save 
the computational costs. Here, these are chosen to be 
free parameters to reproduce a more consistent spectrum 
with the original one. 

Among these seven free parameters, we focus on / cut 
because it depends most strongly on the compactness C 
and the EOS of the NS. Figure [TT1 plots /cut^ici, obtained 
in this fitting procedure, as a function of C. Also the typi- 
cal quasinormal-mode frequencies, /qnm, of the remnant 
BH calculated in Sec. IIVDI are plotted by the two hor- 
izontal lines, which show that the values of / cu tmo for 
models H-Q3M135, HB-Q3M135, and B-Q3M135 agree 
approximately with /qnm and indicates that / cut for 
these models are irrelevant to the tidal disruption. For 
Q = 3, /cut^o depends clearly on the EOS only for 
C < 0.16. This agrees with the result with T = 2 poly- 
tropic EOS (24J. By contrast, / cu twio for Q — 2 depends 
strongly on the NS compactness C irrespective of Mns 
not only for the piecewise polytropic EOS but also for 
r = 2 polytrope [2J]. The solid line in Fig. [11] is the lin- 
ear fitting of ln(/ cut mo) as a function of ln(C) for Q — 2 
and for the piecewise polytrope with r 2 — 3, and denoted 
by a well-approximated relation 

ln(/ cut m ) = (3.87 ± 0.12) InC + (4.03 ± 0.22). (49) 

Thus, fcutirio is approximately proportional to C (for 
Q = 3 and r 2 — 3, /cut^o also appears to be pro- 
portional to C 4 , although the number of data points is 
small and thus this is not conclusive). This is a note- 
worthy point because the power of C is much larger than 
1.5, which is expected from the relation for the mass- 
shedding limit, Eq. (|46[) . Qualitatively, this increase in 
the power is natural because the duration of a NS for the 
survival against tidal disruption after the onset of mass 
shedding is in general longer for a more compact NS due 
to a stronger central condensation of the mass. Equation 
(H3) implies that the ratio / cu t//shod (> 1), where / shcd is 
the frequency of gravitational waves at the onset of mass 
shedding, is larger for the larger values of C. This is the 
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FIG. 11. /cutm-o as a function of C in logarithmic scales. The 
solid line is obtained by a linear fitting of the data for Q — 2 
and F2 = 3. The short-dashed and long-dashed lines show 
approximate frequencies of quasinormal mode of the remnant 
BH for Q — 2 and Q — 3, respectively. 

preferable feature, for an observer of gravitational waves 
from BH-NS binaries who tries to constrain the EOS of 
the NSs, because the dependence of / cu t"io on the EOS 
is enhanced. 

Comparison of the values of /cut^o for models HB- 
Q2M135 (r 2 = 3.0 and C = 0.1718), HBs-Q2M135 
(r 2 = 2.7 and C = 0.1723), and HBss-Q2M135 (r 2 = 2.4 
and C = 0.1741), for which the value of C is approxi- 
mately identical, shows that / cut mo depends also on the 
adiabatic index of EOS in the central region, T 2 . The 
reason for this is that the NSs with smaller values of 
r 2 (but with the same value of C) have more centrally 
condensed density profile as can be seen from the value 
of Pmax in Table HH and hence, are less subject to tidal 
disruption (/ cut mo becomes larger). Quantitatively, the 
value of /cut mo increases by ~ 20%, when the value of 
r 2 is varied from 3 to 2.4. This result suggests that it 
may be possible to constrain not only the compactness of 
a NS but also its density profile and detailed function of 
P(p) for the EOS, if gravitational waves emitted during 
the merger of low-mass BH-NS binaries are detected. 



D. Properties of the disk 

If a NS is tidally disrupted before it is swallowed by 
the companion BH, a disk may be formed around the BH. 
Figure [T2l plots the time evolution of the rest mass of the 
material located outside the apparent horizon M r>rAli 
defined by Eq. (|4"Tj) . This shows that most of the material 
is swallowed by the BH soon after the onset of the merger 
(or tidal disruption) within ~ 1 ms, but 1%-10% of total 
rest mass survives around the BH to be a disk, if the tidal 
disruption occurs (see Table [V] which lists the numerical 
results of M r>rAH at the end of the simulations for all 
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FIG. 12. Evolution of the rest mass of the material located outside the apparent horizon, Af r >, AH , with an appropriate time 
shift; in these plots, the time at the onset of the merger is taken as the time origin. The top-left panel shows the results 
for models with Q — 2 and Mns = 1.35M0 for all the EOSs employed in this paper (we note that the simulation for model 
H-Q2M135 unfortunately terminated in the middle of the accretion process due to the electrical outage at our institute). The 
top-right panel shows the results for selected models with Mns = 1.35Mq but with different values of Q. The bottom-left panel 
shows the results for selected models with Q = 2 but with the different NS mass Mns- The bottom-right panel is the same as 
the bottom- left panel except for the normalization of the mass, with respect to the initial rest mass M, . 



the models). 

To clarify that the disk will survive for a time duration 
longer than the dynamical time scale of the system, we 
estimate an accretion time scale. Figure [T^] shows that 
for t — t ro erger ^ 5 ms, M r>rAa for each model behaves 
approximately as C exp(— t/td) where C is a constant 
and id is the accretion time scale which we determine by 
a least-square fitting of M r>rAIl (t) at t — t morgcI s» 10 
ms. The fourth column of Table [V] lists the numerical 
results. It is found that the accretion time scale is always 
longer than the dynamical time scale of the remnant disk 
~ 10 ms, and hence, we conclude that the BH-NS merger 
always forms a long-lived accretion disk, if the disk is 
formed [57| . 

Figure [T2] plots the values of M r>rAli estimated at 
t — imorgcr ~ 10 ms as a function of the NS compact- 
ness C and clarifies the dependence of the disk mass on 



the EOS. The disk mass for model Bss-Q2M135 is esti- 
mated at the end of the simulations, 4.83 ms, because it 
already became very small at that time and we stopped 
the simulation. (We also note that the result for model 
H-Q2M135 is not included in Fig. [13l because the simu- 
lation for this model unfortunately terminated just after 
the disk formation due to the electrical outage at our in- 
stitute.) This figure summarizes the key features as fol- 
lows: (i) for a given mass ratio and for a given adiabatic 
index of the core, T2, the disk mass decreases monotoni- 
cally with the increase of C for M r>rAH < 0.1M Q ; (ii) for 
a given mass ratio and for a given NS compactness, the 
disk mass increases slightly with the increase of T2 ; and 
(iii) the disk mass is highly sensitive to the mass ratio of 
the binary, Q, for a given mass and EOS of the NS. In 
the following, we observe these features from Fig. [T^] in 
detail. 
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FIG. 14. Relation between disk mass M r>rAH and the max- 
imum density, p m ax, estimated at t — i mC rgcr ~ 10 ms. The 
maximum density oscillates with time even in the quasista- 
tionary phase, and we here plot a value averaged in one oscil- 
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The top left panel of Fig. [T2] plots the disk-mass evo- 
lution for binaries with Q — 2, Mns = 1.35M Q and for 
all the EOSs employed in this paper. For this sample, 
since Mns is identical, and we find that the 
disk mass increases monotonically with C~ l (see Table Q] 
for C of each model); the disk mass is larger for a model 
for which the tidal disruption occurs at a more distant 
orbit (i.e., for smaller value of / cut , cf. Fig. QT]). This is 
quite reasonable because the earlier onset of tidal disrup- 
tion helps more materials to remain outside the IS CO of 
the BH. 

Comparison of the results for models HB-Q2M135 
(r 2 = 3.0 and C = 0.1718), HBs-Q2M135 (r 2 = 2.7 



and C = 0.1723), and HBss-Q2M135 (r 2 = 2.4 and 
C = 0.1741) indicates that the disk mass depends not 
only on the compactness of the NS but also on the adia- 
batic index of the core, r 2 ; a higher value of T 2 is prefer- 
able for forming a massive disk. This dependence on T 2 
is consistent with the result reported in [58] ; the NS with 
a larger value of the adiabatic index is more subject to 
tidal disruption (tidal disruption occurs for more distant 
orbital separation). The physical interpretation for this 
result is that the degree of central mass concentration 
for NSs of larger values of the adiabatic index is weaker, 
helping earlier tidal disruption (in other words, we may 
say that the tidal Love number or deformability is larger 
for the larger value of T 2 ). 

The top right panel of Fig. [12] plots the disk- mass evo- 
lution for the NS with the same mass (Mns = 1.35M Q ) 
but with different mass ratio Q = 2 and Q — 3 and with 
HB and 2H EOSs. This, together with Fig. Q31 shows 
that the disk mass depends strongly on the mass ratio, 
in particular for the soft EOS. The reason for this is sim- 
ply that the NS is less subject to tidal disruption for 
a larger BH mass (i.e., for weaker tidal force near the 
ISCO). The present result suggests that the disk mass is 
much smaller than 0.01M© for BH-NS binaries with the 
typical NS mass of M NS = 1.2-1. 35M© and C > 0.16, if 
the BH is nonspinning and Mbh ^ 4M Q . Only for the 
case C < 0.16, the disk mass may be larger than 0.01M Q 
even with a high-mass BH companion. This conclusion 
is in agreement with the previous studies [22H241 ) . 

The two bottom panels of Fig. [12] compare the disk- 
mass evolution for models 2H-Q2M12 and 2H-Q2M135 
and for models HB-Q2M12 and HB-Q2M135. In the left 
panel we plot the disk mass in units of Mq while the 
bottom right panel plots the disk mass in units of M». 
We note that the NS radius depends weakly on the mass 
for 1.2M < M NS < 1.35M for both EOSs, and also the 
mass ratio Q is identical for these models. Nevertheless, 
the disk mass depends strongly on the NS mass except for 
models with stiff 2H EOS as seen in Table fVl it decreases 
with the increase of A/ns- Thus, not the NS radius i?NS 
but C is the key parameter for determining the disk mass. 

Before closing this section, we summarize several key 
properties of the remnant disk. Figure [14] plots the rela- 
tion between M r>rAH and the maximal rest-mass density 



of the remnant disk estimated at t — i n 



10 



Pmax 

ms. This clearly shows a strong correlation between 
two quantities. The value of M r>rAH increases approxi- 
mately linearly with p max for M r>rAH < 0.1M Q , and for 
M r>rAH > 0.01M Q , p max is larger than 4 x 10 11 g/cm 3 . 
Because the density is high and the temperature should 
be also high enough (~ 10 MeV if viscous effects or 
magnetohydrodynamic effects are taken into account [59h - 
l6lj|). neutrinos will be copiously produced in such a disk 
in reality. Because of the high density and temperature, 
the cross section to the nucleon will be large enough 
(~ 10~ 41 cm 2 ) to trap neutrinos inside the disk of nu- 



cleon number density n n = 
m n is nucleon mass 1.66 x 



= p/m n 
lO" 24 , 



; mm 



cm where 
. Therefore, 
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a neutrino-dominated accretion disk will be always pro- 
duced, if BH-NS binaries result in a system composed of 
the BH and surrounding disk of mass larger than O.OIMq. 

E. Properties of the remnant BH 

Table [V] shows several quantities associated with the 
remnant BH such as the mass and spin, in addition to 
the disk mass. Unlike the disk mass, the mass and spin 
of the remnant BH depend weakly on the EOS of the NS. 
For given values of Q and -Mns, the BH mass tends to 
be slightly smaller for stiffcr EOS, primarily because the 
fraction of the NS mass swallowed by the BH is smaller 
(the disk mass is larger). The spin does not show such 
a clear dependence. The reason is that the spin angular 
momentum of the remnant BH is affected by two com- 
peting processes; one is the orbital angular momentum 
dissipation due to gravitational radiation reaction and 
the other is the distribution of the angular momentum 
to the disk surrounding the BH. The former dissipation 
effect is important for the case in which the NS is com- 
pact and the tidal disruption does not occur as stated 
in Sec. IIVBI By contrast, the latter effect is more im- 
portant for the case in which the NS is less compact and 
the tidal disruption occurs in the relatively early stage of 
the inspiral phase. Although the relation AJ > J r >r A H 
(where J r > r AH denotes the angular momentum of disk) 
always holds for all the models, we may also have the 
relation AE > M r>rAH . As a result, a nondimensional 
spin parameter, which may be approximately estimated 
by 

Spin angular momentum 
(Mass) 2 

^ (Jo — A J — J r>rAH ) , . 

~(M -A£-M r>rAH )2' ^ > 

does not depend simply on the EOS. 

The spin of the remnant BH is primarily determined 
by the mass ratio, Q; a = 0.66 ± 0.03 for Q = 2 and 
a = 0.54 ± 0.02 for Q — 3 (here ± signs do not imply the 
error bars but signify differences due to the EOS). Thus, 
the spin parameter is modified by the EOS only in ±5%. 

From the typical value of the spin parameter a and 
mass of the remnant BH Mbh,^ we estimate typical 
quasinormal-mode frequencies /qnm of the remnant BH 
by the latest fitting formula [65j 

/qnmMbhj « i[1.5251 - 1.1568(1 - a) ' 1292 ]. (51) 

Then, / QNM « 0.083/M BH ,f for Q = 2 and « 
0.076/AfBH,f for Q — 3, respectively. Assuming that 
C e /47r gives an approximate value of AfeH.f as described 
in Sec. IIIIB1 these values are in good agreement with 
the ringdown part of gravitational waves for models Bss- 
Q2M135 and B-Q3M135, for which the disk masses are 
negligibly small, respectively. We note that this estima- 
tion is valid only when the quasinormal modes of the BH 



are excited, and actually the tidal disruption of the NS 
often suppresses the quasinormal-mode excitation as can 
be seen in Figs. [6] and [3 in particular, for the stiff EOS 
such as 2H. 



V. SUMMARY 

We performed numerical simulations for the merger of 
nonspinning BH-NS binaries using an AMR code SACRA 
with eight piecewise polytropic EOSs. In this work, we 
employed the EOSs with two free parameters which de- 
termine the core EOSs. The crust EOS was fixed whereas 
the core EOS was varied for a wide range, to investigate 
the dependence of gravitational waveforms, merger pro- 
cess, and merger remnant on the core EOS. We focused, 
in particular, on the case in which the NS is tidally dis- 
rupted by the companion BH, choosing relatively low val- 
ues of mass ratio as Q — 2 and 3 as well as low masses 
for the NS as Mns = 1-2 and 1.35M Q . By preparing 
the initial condition with a distant orbit and a small ec- 
centricity, we always tracked > 5 quasicircular orbits in 
the inspiral phase and studied the merger phase with a 
realistic setting. We also evolved the merger remnant 
(BH-disk system) until they settled to a quasistationary 
state. 

A wide variety of simulations were systematically 
performed to investigate the dependence of the tidal- 
disruption process and resulting gravitational waveforms 
on the EOS. For the case in which the tidal disruption oc- 
curs before the orbit reaches the ISCO, the gravitational- 
wave amplitude decreases quickly at its onset and the 
emission of ringdown gravitational waves associated with 
the quasinormal mode of the remnant BH is suppressed. 
Only in the BH-NS binaries with low values of mass ra- 
tio (for the nonspinning BH), the tidal effects play an 
important role, and hence, the remarkable dependence 
of the gravitational waveforms on the EOS is found only 
for such cases: With stiffcr EOSs, the radius of the NS 
becomes larger and the tidal effect is more relevant than 
with softer EOSs. For given masses of the BH and NS, 
the tidal disruption occurs in a lower frequency with 
stiffer EOSs than with softer EOSs, and consequently, the 
emission of gravitational waves terminates at a lower fre- 
quency in the inspiral phase. The corresponding Fourier 
spectrum of gravitational waves is characterized by a 
cutoff frequency, / cu t, above which the spectrum am- 
plitude exponentially damps. From the analysis of the 
gravitational-wave spectra, we find that the cutoff fre- 
quency / cu t depends strongly on the mass ratio and the 
compactness C of the NS. For a given small mass ratio 
such as Q = 2, the value of / cu t increases monotonically 
and steeply with C, depending weakly on the adiabatic 
index, T2, of the core EOS. We derive the relation be- 
tween C and / cu t for Q = 2 and r2 = 3 as / cu t oc C 3 9 , in 
which the power index of C is significantly larger than 1.5 
which is expected from the analysis of the mass-shedding 
limit. This implies that the dependence of / cu t on C is 
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TABLE V. Several key quantities for the merger remnants. All the quantities are estimated when we stopped the simulation at 
t = tend- ^merger denotes the time of the merger and the time duration for following the disk evolution, i cnt j — t metgeT , is shown 
in the second column. M r>rAH is the rest mass of the disk surrounding the BH; because the accretion is still ongoing at the 
end of simulations due to the hydrodynamic angular momentum transport process, the values listed give only an approximate 
mass of the long-lived accretion disk (especially for model H-Q2M135; see Sec. IIVD[) . which survives for a time scale longer 
than the dynamical time scale ~ 10 ms. id is the approximate accretion time scale estimated around ~ 10 ms after the merger, 
which we show only for the case M r>TAH > 0.001M©. C e and C p are the circumferential radii of the apparent horizon along 
the equatorial plane and meridional plane, respectively, and C e /47r is the approximate mass of the remnant BH. Mj rr is the 
irreducible mass of the remnant BH. a is the nondimensional spin parameter of the remnant BH estimated from C p /C e . 



Model t on d 


^merger (ills) 


M r>rAH [M e ] 


t d (ms) 


C e /47rM 


Mirr/Mo 


Cp/Ce 


a 


2H-Q2M135 


12.6 


0.097 


30 


0.947 


0.891 


0.913 


0.64 


H-Q2M135 


4.22 


0.070 




0.968 


0.905 


0.905 


0.66 


HB-Q2M135 


11.9 


0.022 


18 


0.978 


0.912 


0.902 


0.67 


HBs-Q2M135 


13.7 


0.015 


15 


0.980 


0.914 


0.902 


0.67 


HBss-Q2M135 


13.5 


0.0093 


12 


0.981 


0.916 


0.903 


0.67 


B-Q2M135 


20.0 


0.0045 


20 


0.980 


0.917 


0.905 


0.66 


Bs-Q2M135 


12.5 


0.0029 


19 


0.979 


0.917 


0.906 


0.66 


Bss-Q2M135 


4.83 


6 x 10 -4 




0.977 


0.917 


0.910 


0.65 


2H-Q3M135 


20.4 


0.044 


19 


0.961 


0.925 


0.944 


0.52 


H-Q3M135 


21.4 


0.0015 


11 


0.984 


0.942 


0.937 


0.55 


HB-Q3M135 


16.5 


2 x 10~ 4 




0.983 


0.942 


0.937 


0.55 


B-Q3M135 


15.1 


< 10~ 5 




0.982 


0.941 


0.939 


0.55 


2H-Q2M12 


12.6 


0.097 


31 


0.939 


0.886 


0.918 


0.62 


H-Q2M12 


11.9 


0.077 


30 


0.959 


0.899 


0.907 


0.66 


HB-Q2M12 


9.72 


0.068 


30 


0.964 


0.902 


0.906 


0.66 


B-Q2M12 


15.4 


0.043 


24 


0.972 


0.908 


0.903 


0.67 


HB-Q3M12 


12.2 


0.011 


15 


0.979 


0.937 


0.937 


0.55 


B-Q3M12 


10.6 


0.0019 


17 


0.982 


0.940 


0.936 


0.56 



stronger than that for /shod, and indicates that the obser- 
vation of / cu t will play a role for constraining the value of 
C. Varying the core EOS also modifies the value of / cut , 
because the central density profile of the NS depends on 
the stiffness of the core EOS and susceptibility to the 
tidal force of its companion BH is modified. For the vari- 
ation from r2 = 3 to 2.4, the value of / cu t is modified 
by ~ 20%. This suggests that the details of the core 
EOS for p > 10 15 g/cm 3 may play an important role for 
determining the gravitational waveform from the BH-NS 
binaries composed of high-mass NSs. 

We also determined the mass of the disk surrounding 
the remnant BH. The disk mass depends strongly on the 
EOS, because the EOS determines the location at which 
the tidal disruption occurs through the compactness C 
of the NS. The disk mass is correlated strongly with the 
NS compactness C, and for Q = 2, it can be > O.OlAi© 
for a wide range of the EOSs and the NS masses Mns- 
However, the disk mass is tiny for Q = 3, unless the EOS 
is extremely stiff like 2H EOS or the NS mass is low. For 
the BH-NS binaries consisting of a nonspinning BH, the 
disk mass can be > 0.01M Q for Q = 3, only for the case 
C<0.16. 

Using the quantities calculated on the apparent hori- 
zon, wc estimated the dimensionless spin parameter of 
the remnant BH. We find that this spin parameter de- 



pends only weakly on the EOS for given masses of the 
BH and NS, unlike the disk mass. The BH spin depends 
primarily on the mass ratio Q and becomes smaller for 
a binary with a larger value of Q: aw 0.66 ± 0.03 for 
Q = 2 and « 0.54 ± 0.02 for Q = 3. 

Finally we list the issues for the future. The two-piece 
EOS employed in this paper is not accurate enough to 
describe high-mass NSs for which the inner core is com- 
posed of a high-density matter with p > 10 15 g/cm 3 . For 
the study of a BH-NS binary composed of a high-mass 
NS with small values of Q (i.e., for a binary in which the 
tidal interaction plays a role), it is necessary to adopt 
piecewise polytrope EOSs with three or four free param- 
eters. It is also necessary to take into account the BH 
spin for a systematic survey of the BH-NS binary merger 
process, because the orbital frequency at the ISCO de- 
pends strongly on the BH spin as well as the mass of the 
BH; e.g., for the ISCO around Kerr BHs, the orbital an- 
gular frequency increases by a factor of 6 3 / 2 if the spin is 
changed from zero to unity. This difference in the ISCO 
will be crucial for determining the criteria for the onset 
of tidal disruption, the mass of the remnant disk, and 
gravitational waveforms. Currently we are working on 
this subject and will report the numerical results in the 
next paper. 
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TABLE VI. Several numerical results for models HB-Q2M135 
and H-Q3M135 with different grid resolutions, N = 50, 42, 
and 36. All the quantities are defined in the body text. In 



this table, we compare the disk mass at t — t n 



10 ms. 



N /cutm M r>rAH [M ](lOms) 


a 


AE/M (%) AJ/M%) 


HB-Q2M135 


50 0.0613 


0.025 


0.67 


1.36 


21.8 


42 0.0621 


0.022 


0.67 


1.34 


21.4 


36 0.0644 


0.022 


0.68 


1.35 


21.4 


H-Q3M135 


50 0.0790 


0.0027 


0.55 


1.39 


21.8 


42 0.0794 


0.0021 


0.56 


1.36 


21.3 


36 0.0788 


0.0022 


0.56 


1.33 


20.7 



0.2 



0.1 
0.08 

E 0.06 
q 0.05 
£ 0.04 

'£ 0.03 
0.02 



0.01 




HB-Q2M135 N=50 
HB-Q2M135 N=42 
HB-Q2M135 N=36 
H-Q3M135 N=50 
H-Q3M135 N=42 
H-Q3M135 N=36 



0.01 



0.02 0.03 0.04 0.050.06 0.08 0.1 
f m n 
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APPENDIX: CONVERGENCE 



FIG. 15. Comparison of gravitational-wave spectra for models 
HB-Q2M135 and H-Q3M135 with different grid resolutions. 
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In this Appendix, we demonstrate that the convergence 
is approximately achieved for the numerical results shown 
in Sec. HVl We here compare numerical results obtained 
with different grid resolutions for models HB-Q2M135 
and H-Q3M135. Table IVTl lists several numerical results. 
This shows that the numerical results depend only weakly 
on the grid resolutions, and thus, we conclude that the 
convergence is approximately achieved in our simulation. 
Most importantly, Fig. [15] shows that the gravitational- 
wave spectra approximately converge and /cut^o shown 
in Table |VT] does not vary by > 5%. Figure [16] plots 
/cut^io for model HB-Q2M135 as a function of the inverse 
of a squared grid resolution 1/N 2 . This figure shows that 
the value of / cut mo converges at better than second order, 
and thus the values of / cu t?7io for = 50 are obtained in 
< 3% error. For model H-Q3M135, the value of fcutm-o 
does not converge systematically and fluctuates with the 
amplitude of ~ 0.5%. This fluctuation may be ascribed 
to the variance associated with the fitting procedure us- 
ing Eq. (|47[) . which involves a number of free parameters. 
We estimate roughly the variance of / cu tTOo at ~ 0.5% 
within 95% accuracy of the fitting with respect to the 



FIG. 16. /ctTOo as a function of the inverse of a squared grid 
resolution 1/N 2 for model HB-Q2M135. 



norm defined by Eq. gSJl for model H-Q3M135. We note 
that the merger time t mergel depends on the grid resolu- 
tion; it is systematically larger for better grid resolutions. 
However, the spectrum near / = / cut depends weakly on 
the grid resolution. AE and A J also approximately con- 
verge. The errors are < 0.1% for AE and < 1% for A J, 
respectively. 

Among many quantities, the disk mass is most sensi- 
tive to the numerical dissipation because the spurious dis- 
sipation of the angular momentum in the disk enhances 
the accretion of the materials surrounding BH and results 
in a lower disk mass. Hence, the values of disk mass de- 
scribed in the body text should be regarded as the lower 
limit of the actual mass of the remnant disk. 

We plot the time evolution of the disk mass for dif- 
ferent grid resolutions in Fig. [17] Roughly speaking, the 
numerical results for the disk mass increase with improv- 
ing the grid resolution, although systematic convergence 
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FIG. 17. Comparison of the disk-mass evolution for models 
HB-Q2M135 and H-Q3M135 with different grid resolutions. 

property is not seen. The reason for this unsystematic 
behavior is likely that the motion of the disk material is 
affected slightly by the atmosphere (in particular for low- 
mass disks of relatively low densities), and thus, conver- 
gence property should not be expected. Assuming most 
conservatively that the convergence is achieved only at 
first order for the results of N — 42 and 50, the error of 
the results for N = 50 may be a factor of 2 for the low 
disk mass case M r>rAU < 0.01Af Q . However, we expect 
that systematic quantitative relations between the disk 
mass and the compactness of the NS, and between the 
disk mass and the maximum density shown in Figs. 1131 
and Q3] are not drastically changed. 



VI. ERRATUM 

(The original version of this article was submitted on 
Aug 9, 2010. This erratum is added on Aug 5, 2011.) 

In previous sections, we performed numerical simula- 
tions for the merger of a nonspinning black hole (BH) and 
a neutron star (NS), and explored gravitational waves 
emitted and the final outcome formed after the merger. 
We recently noticed that we systematically underesti- 
mated disk masses in previous sections. The reason is 
that we evolved hydrodynamic variables and estimated 
disk masses only in domains of the size ~ 200 3 km 3 , al- 
though Einstein's field equation was solved in domains 
of the size ~ 800 3 km 3 . A small domain size for hy- 
drodynamics is insufficient for the estimation of the disk 
mass, because, if tidal disruption occurs at a distant or- 
bit, especially for the case in which the NS radius is large 
(~ 15 km), tidal disrupted material extends far away 
from the central region. For this reason, we performed 
again the same simulations as in previous sections, en- 
larging the computational domain of hydrodynamics. To 
estimate disk mass more accurately, in addition, we en- 
larged the size of the computational domain to the size 



TABLE VII. Setup of the grid structure for the computation 
with our AMR algorithm. Ax = h 7 = L/(2 7 N) is the grid 
spacing at the finest-resolution domain with L being the loca- 
tion of the outer boundaries for each axis. i?diam/A:r denotes 
the grid number assigned inside the semimajor diameter of 
the NS. Ao is the gravitational wavelength of the initial con- 
figuration. 



ivioQei 


L\X j IVlQ 


-Kdiam/ L\X 


Xjf AO 


2H-Q2M135 


0.0471 


90.8 


2.377 


H-Q2M135 


0.0377 


86.2 


2.130 


HB-Q2M135 


0.0347 


87.0 


1.963 


HBs-Q2M135 


0.0353 


85.2 


1.996 


HBss-Q2M135 


0.0353 


84.0 


1.996 


B-Q2M135 


0.0330 


85.1 


1.863 


Bs-Q2M135 


0.0324 


84.4 


1.830 


Bss-Q2M135 


0.0270 


95.4 


1.650 


2H-Q3M135 


0.0353 


89.0 


1.996 


H-Q3M135 


0.0282 


84.7 


1.711 


HB-Q3M135 


0.0269 


82.7 


1.631 


B-Q3M135 


0.0247 


83.8 


1.497 


2H-Q2M12 


0.0565 


86.9 


2.510 


H-Q2M12 


0.0453 


83.1 


2.563 


HB-Q2M12 


0.0420 


83.6 


2.377 


B-Q2M12 


0.0392 


83.4 


2.218 


HB-Q3M12 


0.0306 


84.6 


1.713 


B-Q3M12 


0.0278 


86.9 


1.572 



1500 3 -2000 3 km 3 . This is done by increasing a coarse do- 
main by one more level in the adaptive mesh refinement 
algorithm (AMR). Specifically, the number of coarser do- 
mains is increased from three to four. Table I VIII (new 
version of Table III) summarizes the parameters of the 
new grid structure. In these simulations, the total rest 
mass of the atmosphere is always less than 1O _4 M0. Re- 
sults for gravitational waves do not change within the 
level of numerical accuracy in our simulations. 

Table IVTnl (corrected version of Table V) lists cor- 
rected values for quantities associated with the merger 
remnants. We estimated all the values at the end of 
the simulations in previous sections. In this section, we 
present the values evaluated at 10 ms after the merger 
to perform more systematic comparisons. Quantities as- 
sociated with the remnant BH do not change appreciably. 
Taking into account the change in the time at which the 
disk mass is estimated, the mass of the remnant disk be- 
comes larger by a factor of ~ 2-3 for Q = 2, and by a 
factor of ~ 5 for Q — 3. Figure [18] (corrected version of 
Fig. 12) plots the time evolution of M r>rAii . Although 
qualitative behavior is not altered, M r>rAH is system- 
atically larger for the new computations. In particular, 
the sudden decrease of M r>rAH at ~ 3-5 ms after the 
merger seen in Fig. 12 now disappears. Approximate ac- 
cretion time scale td becomes longer by a factor of < 2 
for many cases. Figure Q1J] (corrected version of Fig. 13) 
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plots the values of M r>rAH at « 10 ms after the merger 
as a function of C. Although we again see the systematic 
increase of M r>rAH , the conclusion that the disk mass is 
much smaller than O.OlAf© for BH-NS binaries with the 
typical NS mass of M NS = 1.2-1.35M and C > 0.16 
does not change. Figure [20] (corrected version of Fig. 14) 
plots the relation between M r>rAM and the maximum rest 
mass density /o ma x of the remnant disk. Approximately 



speaking, the relations between them are not changed 
qualitatively and quantitatively. 

Table IIXI (corrected version of Table VI) lists several 
numerical results for the merger remnants, and Fig. [5T] 
(corrected version of Fig. 17) plots the time evolution of 
M r>rAH for different grid resolutions. The convergence 
of the remnant disk mass becomes slightly better than 
that in previous sections. 
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TABLE VIII. Several key quantities for the merger remnants. All the quantities are estimated at t — tmerger ~ 10 ms, where 
^merger denotes the time of the merger. M r>rAH is the rest mass of the disk surrounding the BH; because the accretion is still 
ongoing at the end of simulations due to the hydrodynamic angular momentum transport process, the values listed give only 
an approximate mass of the long-lived accretion disk, which survives for a time scale longer than the dynamical time scale ~ 10 
ms. td is the approximate accretion time scale estimated around w 10 ms after the merger, which we show only for the case 
M r>TAH > O.OO1M0. C e and C p are the circumferential radii of the apparent horizon along the equatorial plane and meridional 
plane, respectively, and C e /47r is the approximate mass of the remnant BH. Mi lr is the irreducible mass of the remnant BH. a 
is the nondimensional spin parameter of the remnant BH estimated from C p /C e . 
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TABLE IX. The disk masses at t — i mor gcr ~ 10 ms and nondi- 
mensional spin parameters of the remnant BHs for models 
HB-Q2M135 and H-Q3M135 with different grid resolutions, 
N = 50, 42, and 36. 
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FIG. 18. Evolution ol the rest mass of the material located outside the apparent horizon, M r >,- AH , with an appropriate time 
shift; in these plots, the time at the onset of the merger is taken as the time origin. The top-left panel shows the results for 
models with Q — 2 and A^ns = 1.35M© for all the EOSs employed in this paper The top-right panel shows the results for 
selected models with Mns = 1.35Mq but with different values of Q. The bottom-left panel shows the results for selected models 
with Q — 2 but with the different NS mass Mns- The bottom-right panel is the same as the bottom- left panel except for the 
normalization of the mass, with respect to the initial rest mass M». 
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FIG. 19. Disk mass M r >r Ali at t 
function of the NS compactness C. 
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FIG. 20. Relation between disk mass Mi 
imum density, pmax, estimated at t — t„ 
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and the max- 
s 10 ms. The 



maximum density oscillates with time even in the quasista- 
tionary phase, and we here plot a value averaged in one oscil- 
lation period. 
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FIG. 21. Comparison of the disk-mass evolution for models 
HB-Q2M135 and H-Q3M135 with different grid resolutions. 



